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ABSTRACT 

We investigate the effect of shear viscosity, is, and Ohmic resistivity, i] on 
the magnetorotational instability (MRI) in vertically stratified accretion disks 
through a series of local simulations with the Athena code. First, we use a se- 
ries of unstratified simulations to calibrate physical dissipation as a function of 
resolution and background field strength; the effect of the magnetic Prandtl num- 
ber, P m = u/i], on the turbulence is captured by ~ 32 grid zones per disk scale 
height, H . In agreement with previous results, our stratified disk calculations are 
characterized by a subthermal, predominately toroidal magnetic field that pro- 
duces MRI-driven turbulence for \z\ < 2H. Above \z\ ~ 2H, magnetic pressure 
dominates and the field is buoyantly unstable. Large scale radial and toroidal 
fields are also generated near the mid-plane and subsequently rise through the 
disk. The polarity of this mean field switches on a roughly 10 orbit period in a 
process that is well-modeled by an a-VL dynamo. Turbulent stress increases with 
P m but with a shallower dependence compared to unstratified simulations. For 
sufficiently large resistivity, rj ~ c s H / '1000 where c s is the sound speed, MRI tur- 
bulence within 2H of the mid-plane undergoes periods of resistive decay followed 
by regrowth. This regrowth is caused by amplification of toroidal field via the 
dynamo. This process results in large amplitude variability in the stress on 10 to 
100 orbital timescales, which may have relevance for partially ionized disks that 
are observed to have high and low accretion states. 
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Introduction 



Disk accretion is a fundamental astrophysical process, responsible for such diverse phe- 
nomena as the immensely luminous, distant quasars and the formation and evolution of 
protostellar systems. This process requires removing angular momentum from the orbiting 
gas, and from the beginnings of disk theory it has been clear that microphysical viscosity is 
orders of magnitude too small to account for observed accretion rates. It is now understood 
that accretion is driven by magnetohydrodynamic (MHD) turbulent stresses arising from a 



robust and powerful instability known as the magnetorotational instability (MRI) (Balbus 



& Hawley 1991, 1998). 



Analytic investigations of the MRI have proven to be very insightful (e.g., Balbus & 
Hawley|199l| |Goodman fc Xu|1994[ ), but they can offer only limited guidance to the behavior 
of the MRI in the fully turbulent saturated state. Numerical simulations have become the 
essential tool for investigating how the MRI operates in accretion disks, and in particular, 
local shearing box simulations (see |Hawley et aL|l995[ ) have proven to be especially useful in 
such studies. The shearing box system solves the equations of MHD in a local, co-rotating 
patch of an accretion disk. The size of this patch is assumed to be small compared to its 
radial distance from the central star, allowing one to use Cartesian geometry while retaining 
the essential dynamics of differential rotation. Shearing box simulations have led to an 
improved understanding of MRI turbulence while addressing some basic questions about 
accretion disks. 

One such question is, what sets the saturation level of MRI-driven turbulence and the 
angular momentum transport rate? This is usually phrased as "What is a?" following 



Shakura & Syunyaev (1973) who made the ansatz that the rcf) component of the turbulent 

aP. Previous shearing box simulations have 



stress, T r( £, is proportional to the pressure, T rq 
found substantial evidence that MRI-driven stress does not behave in a manner consistent 



with how a is often applied in disk theory. The early studies of Hawley et al. (1995) showed 



that the stress is proportional to the magnetic pressure, but the magnetic pressure is not itself 
directly determined by the gas pressure, a result that holds over a wide range of shearing 



box simulations (Blackman et al. 2008). Sano et al. (2004) performed an extensive survey 



and observed, at best, only a very weak gas pressure dependence for turbulent stress. More 



recently, vertically stratified local simulations of radiation-dominated disks (Hirose et al. 



2009) have shown that while there is a correlation between the MRI stress and total (gas 
plus radiation) pressure, it has the opposite causal relationship from that normally assumed 
in the a model: the stress determines the pressure, not the other way around. Increased 
stress can lead to increased pressure through turbulent heating, but an increase in pressure 
does not feed back into the stress. 
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What then does determine the saturation level of the MRI? We know that viscosity and 
resistivity can significantly affect the properties of the linear MRI by reducing growth rates 
and altering the stability limits. Recently, the influence of the viscous and Ohmic dissipation 
on MRI-induced turbulence has become a focus of shearing box simulations. Previous inves- 
tigations into the effect of Ohmic resistivity on the saturated state ( Hawley et al.|1996 Sano 
et al.|1998t [Fleming et al.|2000t[Saiio fc Inutsuka(200T| |Ziegler fc Rudiger|200T| |Sano fc Stone 



2002b) found that angular momentum transport decreases as the resistivity is increased, but 



it was not until the very recent work of Fromang & Papaloizou (2007), Pessah et al. (2007) 



and Fromang et al. (2007) that the influence of physical dissipation in numerical simulations 



was fully appreciated. Fromang & Papaloizou (2007) and Pessah et al. (2007) found that 



for shearing box simulations without a net magnetic flux and without explicit dissipation 
terms, the turbulent saturation level decreases with increasing numerical resolution without 
any sign of convergence. This surprising result was subsequently confirmed with different 
numerical codes 



(e.g. 



Simon et al. 


2009; 


Guan et al. 


2009) 



et al. (2007) demonstrated that the saturation level of the MRI in shearing boxes depends 
strongly on the magnetic Prandtl number, i.e., the ratio of viscosity to resistivity (P m = vjr\). 
Specifically, for local simulations without a net magnetic flux, the turbulent stresses increase 
nearly linearly with P m and for P m < 1, the turbulence decays. Fromang (2010) subsequently 
showed that at a fixed P m , the presence of even a small viscosity and resistivity is sufficient to 
provide convergence in the zero-net field case. The increase in stress with P m is also present 



for net vertical fields (Lesur & Longaretti 2007) and net toroidal fields (Simon & Hawley 



2009); for these field geometries, the P m dependence is weaker and turbulence is maintained 
even for P m < 1. For the net toroidal field model, only a sufficiently high resistivity could 



actually kill the turbulence completely (Simon & Hawley 2009). 



Most of these viscosity and resistivity studies used the unstratified shearing box model 
where vertical gravity is ignored and periodic boundary conditions are assumed for the 
vertical direction. Shearing box models that include vertical stratification are more realistic 
however, and interesting new effects arise when vertical gravity is included. For example, in 
most such simulations, net toroidal field is generated near the mid-plane via MRI turbulence, 
buoyantly rises upwards, and is replaced with a field of the opposite sign in the mid-plane 
region. This behavior happens on a timescale of ~ 10 orbits and appears to be indicative 



of an MHD dynamo operating within the disk (e.g., Brandenburg et al. 1995 Stone et al. 
19961 IHirose et all|2006l |Guan fc Gammie||20T0t |Shi et al.||20T0| |Gressel|[20l"0l |Davis et al. 



2010). Furthermore, the vertical structure of the disk consists of MRI-turbulent gas that is 
marginally stable to buoyancy within \z\ ~ 2H, whereas outside of this region, the gas is 
magnetically dominated, significantly less turbulent, and buoyantly unstable (e.g., Guan & 



Gammie 


2010 


Shi et al. 


2010 
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The effects of resistivity on the MRI in vertically stratified shearing boxes has been 



studied (e.g., Miller & Stone 2000), but many of these calculations are designed with pro- 



tostellar systems in mind and employ a very large resistivity to completely quench the MRI 



and create a "dead zone" near the midplane (e.g., Gammie 1996 Fleming & Stone 2003 



Fromang fc Papaloizou||2006t |Uishi et al.||2007| |Turner fc Sanol|2008| |Ilgner fc NeTson][2008 



Oishi & Low 2009 Turner et al. 2010). The effects of smaller resistivities and of v and P n 



on vertically stratified turbulence has barely been examined. Recently, however, Davis et al. 



(2010) investigated how the results of Fromang & Papaloizou (2007) and Fromang et al. 



(2007) might change in the presence of vertical stratification. Davis et al. (2010) consid- 



ered the zero-net magnetic flux case and employed vertically periodic boundary conditions 
to ensure that zero net flux was maintained throughout the simulation. They found that 
without physical dissipation, the volume-averaged stress level reaches a constant value as 
numerical resolution is increased, in contrast to unstratified simulations where the stress 
declines. They further examined three v and rj values that lead to a decay of turbulence 
in unstratified boxes. With vertical gravity, however, these v and i] values lead to large, 
long timescale fluctuations in the volume averaged stress level; the turbulence saturates at a 
relatively high level for ~ 100 orbits, then decreases to a lower saturation level for another 
~ 100 orbits, and then increases again. 

The primary goal of this work is to extend these first results and investigate in more 
detail how v and rj affect MRI turbulence in vertically stratified shearing boxes. In partic- 



ular, we consider the origin of the long-term fluctuations seen in Davis et al. (2010), and 



whether that effect might be relevant to real accretion disks. Since stratification can al- 
ter the saturation levels of MRI-induced turbulence, we also investigate whether increasing 
P m still leads to increased stress levels and how P m affects the dynamo process previously 
observed in stratified simulations. Finally, these simulations will also serve as an essential 
starting point for future studies that include more realistic physics, such as temperature- 
and density- dependent v and rj. 

The structure of this paper is as follows. In § [2j we describe our evolution equations, 
parameters, and initial conditions. In § [3j we present a series of unstratified shearing box 
simulations to calibrate the effects of physical dissipation and serve as controls for the ver- 
tically stratified shearing boxes with constant v and rj. We discuss our vertically stratified 
simulations in § |4j which are the primary focus of this paper. The first set of these sim- 
ulations contain no physical dissipation, and we carry out several analyses to improve our 
understanding of vertically stratified MRI turbulence. The second set of simulations then 
includes physical dissipation to study the P m effect. We wrap up with a discussion and our 
general conclusions in § |5j We also present a detailed description of our numerical algorithm 
in an Appendix. 
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2. Method 

In this study, we use Athena, a second-order accurate Godunov flux-conservative code 
for solving the equations of MHD. Athena uses the dimensionally unsplit corner transport 



upwind (CTU) method of Colella (1990) coupled with the third-order in space piecewise 



parabolic method (PPM) of Colella & Woodward (1984) and a constrained transport (CT; 
Evans fc Hawley||1988 ) algorithm for preserving the V B = constraint. We use the HLLD 



Riemann solver to calculate the numerical fluxes (Miyoshi & Kusano 2005 Mignone 2007). 



A detailed description of the Athena algorithm and the results of various test problems are 
given in 



Gardiner & Stone (2005), Gardiner & Stone (2008), and Stone et al. (2008). 



Our simulations utilize the shearing box approximation, a model for a local co-rotating 
disk patch whose size is small compared to the radial distance from the central object, R . 
We construct a local Cartesian frame, x = (R — R ), y = R <f), and z, co-rotating with an 
angular velocity Q corresponding to the orbital frequency at R Q , the center of the box. In 



this frame, the equations of motion become (Hawley et al. 1995): 



dp 
di 



+ V • (pv) = 0, 



dpv 



+ V • (pvv - BB) + V [P 



dB 

~dt 



-B' 



2qptt 2 x - ptt 2 z - 2fL x pv + V • T, 



V x (v x B) = -V x (77V x B) 



(1) 
(2) 

(3) 



where p is the mass density, pv is the momentum density, B is the magnetic field, P is the 
gas pressure, and q is the shear parameter, defined as q = —d\nQ/d\iaR. We use q = 3/2, 
appropriate for a Keplerian disk. We assume an isothermal equation of state P = pc 2 , 
where c s is the isothermal sound speed. From left to right, the source terms in equation ^ 
correspond to radial tidal forces (gravity and centrifugal), vertical gravity, the Coriolis force, 
and the divergence of the viscous stress tensor, T, defined as 



T 



pv 



( dvi dvj 



+ 



V dxj dx 



-SijV ■ v 



(4) 



where the indices refer to the spatial components (Landau & Lifshitz 1959), and v is the 



shear viscosity. We neglect bulk viscosity. The source term in equation ^ is the effect of 
Ohmic resistivity, i], on the magnetic field evolution. Note that our system of units has the 
magnetic permeability p = 1. Details about the numerical integration of these equations are 
presented in the Appendix. 
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For unstratified shearing box simulations, the boundary conditions are periodic for y and 
z, and shearing-periodic for x. In stratified simulations, the periodic z boundary conditions 
are replaced with outflow boundary conditions. The specifics are described in the Appendix. 



2.1. Dissipation Parameters 

In all of our simulations, v and r\ are parameterized in terms of Reynolds numbers. 
At a scale height, H, away from the radial center of the shearing box, the fluid velocity is 
1^1 ~ qHQ ~ c s . Thus, the sound speed is a representative velocity for the fluid, and we 
define the Reynolds numbers as 

Re = (5) 
Similarly we define the magnetic Reynolds number, 

Rm ee (6) 
V 

and their ratio, the magnetic Prandtl number, 

P. - - = (7) 
T) Re 

Because resistivity affects the MRI directly, another useful dimensionless quantity is the 
Elsasser number 



A can be computed on a zone-by-zone basis, but it is also helpful to calculate an average A 
for a simulation. Thus, we can write the characteristic Alfven speed in direction % via the 
averaged magnetic field in that direction, 



where the angled brackets denote a volume average, and the subscript i = (x, y, z) depending 
on the direction of interest. In the above definition of A, all three components are used in 
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calculating the Alfven speed. As the Elsasser number approaches unity, resistivity becomes 



more dominant, stabilizing the MRI. Simon & Hawley (2009) found that turbulence decayed 



in net toroidal field, unstratified shearing boxes when the Elsasser number was < 100. 

Although we focus on physical dissipation, there is numerical dissipation as well. One 
way to characterize the effects of finite resolution is to compare the characteristic MRI 
wavelength, Amri = 2tiva/Q to the grid zone size. Noble et al. (2010) defined this "quality 
factor" Q as 



Qi = 



A 



MRI,i 
Ax,; 



27lVAi 

QAxj 



(10) 



Again, vm is defined via equation (|9]). For Qi < 6, the growth of the MRI can be reduced 
( Sano et al.|2004 ) and the MRI is under-resolved, although this number has some uncertainty 
and should be taken only as an estimate. 

One can also calculate a volume-averaged Alfven speed rather than an Alfven speed 
determined by the averaged field, i.e., 



VAi 



\Bi[ 
y/P 



'ir. 



For convenience, most of our analysis uses equation (|9]), as the volume-averaged history data 
and the one-dimensional, horizontally averaged quantities are routinely computed at high 
time resolution for use in creating, e.g., space-time plots. As a check, we have calculated the 



volume- averaged Alfven speed via equation ( 11 ) for several hundreds of orbits in a few of our 



vertically stratified simulations (the volume average is done within 2H of the mid-plane). We 
found that Va from equation (11) is roughly 1-1.5 times larger than that from equation Q. 



2.2. Parameters and Initial Conditions 

We have run shearing box simulations with different field geometries, dissipation values 
and resolutions, both with and without vertical gravity. Here, we describe the parameters 
and initial conditions used in these two types of simulations. The results from the unstratified 
simulations are presented in § [3] and the results from the stratified simulations are in § |4} 
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2.2.1. Sine Z Unstratified Simulations 



The first set of simulations are unstratified, zero net magnetic flux shearing boxes. This 



is the same type of problem as investigated in Fromang et al. (2007) and Simon & Hawley 



(2009). The initial magnetic field is 



B = ^2P /P sm(2nx/L x ) 



(12) 



with (3 = 400. The isothermal sound speed is c s = 0.001, corresponding to an initial gas 
pressure P Q = 10~ 6 with initial density p Q = 1. The orbital velocity of the local domain is 
Q = 0.001. For these simulations, we define the scale height to be 



H 



n 



(13) 



which is a slightly different definition than that in the stratified box (by a factor of y/2). 
The size of the box is L x = 1H, L y = 4H, and L z = 1H. We run several resolutions in order 
to study convergence, and at each resolution, we study four different cases of Rm and P m 
values. See Table [T] for the parameters of these runs. The labeling scheme of the runs refers 
to resolution, field geometry, and dissipation values; e.g.,16SZRml2800Pml6 corresponds to 
16 zones per H, "SZ" for the "Sine Z" geometry, and Rm = 12800, P m = 16. Since, as we 
will see, Rm is the critical parameter in many of our simulations, we choose to include it 



as part of the run label, differing from the convention of related works (e.g., Fromang et al. 



2007 Simon & Hawley 2009) 



The MRI is seeded with random perturbations to the density and the velocity compo- 
nents introduced at the grid scale. The amplitude of the density perturbations is Sp = 0.01 
and the amplitude of the velocity perturbations is (l/5)5pc s for each component (a differ- 
ent perturbation is applied for each component). We do not employ orbital advection in 
these simulations (see Appendix). All simulations are run to 400 orbits, except for the runs 
in which the turbulence decays and also 32SZRml2800Pml6 and 32SZRel2500Pm4, which 
were run to 289 orbits and 246 orbits, respectively. We also include the results from previous, 



higher resolution simulations (Simon & Hawley 2009) in our analysis. 



2.2.2. Flux Tube Unstratified Simulations 



The second set of unstratified simulations contain a net toroidal field and are initialized 
with the twisted azimuthal flux tube of Hirose et al. (2006), with minor modifications to 
the dimensions and /3 values. The initial conditions are designed to match those used in the 
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vertically stratified simulations of § 2.2.3 so we define H to be that of a stratified, isothermal 
disk, 

(H) 

The isothermal sound speed, c s = 7.07 x 10~ 4 , corresponding to an initial value for the gas 
pressure of P Q = 5 x 1CT 7 . With Q = 0.001, the value for the scale height is H = 1. As in the 
SZ simulations, random perturbations are added to the density and velocity components. 

The initial toroidal field, B y , is given by 



" ifB| + B« 2 = 

and we have run simulations with a toroidal field (3 of (3 y = 100, 1000, and 10000. The 
poloidal field components, B x and B z , are calculated from the y component of the vector 
potential, 

Av = /-^£ [! + «»(¥)] ifr<f 
I ifr>f 
where r = a/x 2 + ^ 2 and /? p = 1600 is the poloidal field (3 value. 



The labeling convention is the same as in § 2.2.1, but with "FT" for "Flux Tube" (see 
Table Q. We run simulations both with and without physical dissipation. Those simulations 
with no physical dissipation are labeled with "Num" for "Numerical dissipation", and for 
the y = 1000 (10000) cases, we append /31000 (£10000) on the end of the run label. The 
domain size is L x = 2H, L y = 4H, and L z = 1H, and the resolution is 32 zones per H. 
Orbital advection is employed in these calculations (see Appendix). 32FTNum was run to 
150 orbits, and both 32FTNum/31000 and 32FTNum/310000 were run to 110 orbits. 

The runs with physical dissipation are initiated from the turbulent state (at t = 100 
orbits) of the corresponding (3 y value run with only numerical dissipation. These simulations 
were all run out to 220 orbits. 



2.2.3. Vertically Stratified Simulations 



In vertically stratified shearing boxes, the initial density corresponds to an isothermal 
hydrostatic equilibrium 
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Table 1. Unstratified Simulations 



Label 


Re 


Rm 


Pxn 


Resolution 
(zones per H) 


a 


Description 


1 fi^7Rml 9£0nPm1 (\ 
±UQZjJTvni J_ZOUUJr III J. u 


son 
ouu 


i 9«oo 


1 

1U 


I (-, 


01 1 


zero net tlnx 


^9S7Rm1 980nPm1 fi 

UZrOZjJ_\,lll A-ZiO\J\JL 111J.U 


800 


12800 


16 


32 


0.033 




U^oZjaxIII-IZoUU-T Iill u 


soo 
ouu 


1 9800 


I P, 

1U 


f, 1 


049 


zero net flux 


128SZRml2800Pml6 a 


800 


12800 


16 


128 


0.046 


zero net flux 


16SZRml2500Pm4 


3125 


12500 


4 


16 


0.0043 


zero net flux 


32SZRml2500Pm4 


3125 


12500 


4 


32 


0.013 


zero net flux 


64SZRml2500Pm4 


3125 


12500 


4 


64 


0.015 


zero net flux 


128SZRml2500Pm4 a 


3125 


12500 


4 


128 


0.013 


zero net flux 


16SZRm6250Pml 


6250 


6250 


1 


16 


decay 


zero net flux 


32SZRm6250Pml 


6250 


6250 


1 


32 


decay 


zero net flux 


64SZRm6250Pml 


6250 


6250 


1 


64 


decay 


zero net flux 


16SZRm25600Pm2 


12800 


25600 


2 


16 


0.0077 


zero net flux 


32SZRm25600Pm2 


12800 


25600 


2 


32 


0.010 


zero net flux 


64SZRm25600Pm2 


12800 


25600 


2 


64 


0.0078 


zero net flux 


32FTNum 








32 


0.021 b 


flux tube, num. dissipation, (3 = 100 


32FTNum/31000 








32 


0.020 c 


flux tube, num. dissipation, = 1000 


32FTNurm310000 








32 


0.018 c 


flux tube, num. dissipation, f3 = 10000 


32FTRm800Pm0.5 


1600 


800 


0.5 


32 


decay 


restarted from 32FTNum 


32FTRm3200Pm0.5 


6400 


3200 


0.5 


32 


0.0094 


restarted from 32FTNum 


32FTRm3200Pm2 


1600 


3200 


2 


32 


0.018 


restarted from 32FTNum 


32FTRm3200Pm4 


800 


3200 


4 


32 


0.028 


restarted from 32FTNum 


32FTRm6400Pm4 


1600 


6400 


4 


32 


0.029 


restarted from 32FTNum 


32FTRm6400Pm8 


800 


6400 


8 


32 


0.041 


restarted from 32FTNum 


32FTRml600Pml/31000 


1600 


1600 


1 


32 


decay 


restarted from 32FTNum/31000 


32FTRm3200Pm2/3 1000 


1600 


3200 


2 


32 


0.015 d 


restarted from 32FTNum/31000 


32FTRm3200Pm2/310000 


1600 


3200 


2 


32 


decay 


restarted from 32FTNum/3 10000 



a These runs were taken from I Simon fc Hawley| flo&fy 
b Time averaged from orbit 20 to 150 
c Time averaged from orbit 20 to 110 
d Time averaged from orbit 120 to 400 
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z 2 \ 



p(x, y, z) = p exp y-jpj > ( 17 ) 

where p a = 1 is the mid-plane density, and H is the scale height in the disk, as defined in 
2.2.2 A density floor of 10~ 4 is applied to the physical domain as too small a density leads 



to a large Alfven speed and a very small timestep. Furthermore, numerical errors in energy 
make it difficult to evolve regions of very small plasma (3. All other parameters and initial 



conditions are identical to the corresponding FT runs of § 2.2.2 All the stratified runs use 
(3 y = 100 for the initial toroidal field strength. 

The domain size for all vertically stratified simulations is L x = 2H, L y = 4H, and 
L z = 8H, and we have run most simulations at 32 grid zones per H, with a few at 64 zones 
per H. Orbital advection is employed in all of these runs (see Appendix). The runs are listed 
in Table [2] The labeling scheme is the same as for the unstratified simulations, but we omit 
the "FT" label as all the stratified runs use the flux tube initial conditions. We study the 
effect of physical dissipation by restarting the equivalent resolution, numerical-dissipation- 
only run at t = 100 orbits. 



3. Calibration of Physical Dissipation in Unstratified Disks 

While our focus is on the effects of physical dissipation on MRI-driven turbulence in 
stratified disks, we have carried out a series of unstratified shearing box simulations to 
address several points. First, what resolution is needed to capture the influence of physical 
dissipation, and how does a given resolution influence the effects due to dissipation terms? 
Second, what is the effect of physical dissipation on different initial background toroidal field 
strengths? As we will see, this last question has direct relevance to our stratified simulations, 
in which the net toroidal field within a given region changes in time via shear and buoyancy. 

Table [T] is a list of the unstratified shearing box simulations. The column labeled "Res- 
olution" lists N x , the number of zones in one H. The column labeled a gives the averaged 
stress normalized by the averaged gas pressure, 

_ ((pvjvy - B x B y )) 

— w>) — ' {) 



where the double brackets denote a time and volume average, the volume average is calculated 
over the entire simulation domain, and the time average is calculated from orbit 20 to the 
end of the run. Since the gas is isothermal, (P) = {p)c^. 
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3.1. Resolving Physical Dissipation 



We begin with a resolution study of unstratified simulations with physical dissipation. 
Figure [T] shows a as a function of resolution for three different P m values in the SZ simulations. 
The N x = 128 data listed in the table and plotted in the figure were taken from |Simon 



Hawley (2009). Note that those runs were similar to, but not the same as those done for 
this paper. Here, the grid zone size is equal in all directions, Ax = Ay 

Az 



Simon & Hawley (2009) runs had Ax 



Az, whereas the 
0.4 Ay. Furthermore, the calculations done 
in Simon & Hawley (2009) were performed with the Roe method for the Riemann solver, in 
contrast to the HLLD solver used here. 

For the SZ simulations that have sustained turbulence, a appears to be converging with 
resolution. More specifically, by 32 grid zones per H, a is within a factor of ~ 1.4 of the 
corresponding value at 128 zones per H. In contrast, a in the P m = 16 and P m = 4 runs 
increases by about a factor of 3 going from 16 to 32 zones per H. The P m = 2 case shows 
a much smaller change and the turbulence dies out for all resolutions with Re = 6250 and 

Prr 



1 in agreement with the higher resolution simulations of Fromang et al. (2007). 



Next, we explore the influence of physical dissipation terms on the FT initial field 
configuration, using 32 zones per H. The time history of the volume-averaged stress for 
these runs is displayed in Fig. [2] There is a clear dependence on the dissipation parameters 
and on P m in particular. For large enough resistivity (i.e., low Rm), the turbulence decays; 
the critical Rm value is ~ 1000, in agreement with the higher resolution simulations of [Simon 



& Hawley (2009). In Fig. p3l we plot the time-averaged a values, averaged from orbit 120 to 



the end of the simulation, versus P m for these runs (asterisks). Note that 32FTRm3200Pm4 
and 32FTRm6400Pm4 have different Rm values, but the same P m and nearly the same 
saturation level. We also plot a from the higher resolution simulations of Simon fc Hawley 



(2009) (see their Table 3). The dashed lines are linear fits to the data in log-log space. 
Assuming S in a oc P^ we find S = 0.54 for 32 grid zones per H, and 5 = 0.33 for 128 
grid zones per H; the P m dependence is steeper at the lower resolution. Furthermore, all of 
the a values for the higher resolution simulations are larger than the corresponding lower 
resolution simulations. 

Both the zero net flux and net toroidal flux results suggest that moderate resolutions, 
i.e., 32 grid zones per H, may be sufficient to capture the general effects of changing v and 
77, at least for the range of Re, Rm, and P m values considered here. This is not to say that 
everything is converged at this resolution. Indeed, Fig. [3] shows a noticeable resolution effect. 
Full convergence likely requires a sufficiently high resolution to ensure that the effective 
numerical dissipation scale is below the viscous and resistive dissipation scales (see e.g., 
Fromang et al.|20"07 Simon et al. 2009 Simon fc Hawley|20 09), but the general dependence of 
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a on dissipation parameters appears to be captured even at 32 zones per H. This conclusion 
is also supported by the recent results of Flaig et al. (2010), which show that ~ 30 zones per 
H in the vertical direction may be sufficient to characterize the turbulence in their stratified 
shearing box simulations. This is an important point, as available computational resources 
limit us to this resolution in carrying out the comprehensive study of physical dissipation 
effects presented in this work. 



3.2. Dissipation and Initial Field Strength 

Because the fastest growing MRI wavelength is proportional to va, the background field 
strength can play a significant role in the outcome of a shearing box simulation. For example, 



Longaretti & Lesur (2010) found that the dependence of angular momentum transport on 
Rm and P m becomes steeper for weaker background (vertical) fields in unstratified shearing 
boxes. 

Here we consider the effects of different toroidal field strengths on the FT initial con- 
dition in the unstratified shearing box. We consider (3 y = 1000 and (3 y = 10000 runs both 
with and without physical dissipation (see Table [l] and §[2]). The time history of the volume- 
averaged stresses for these simulations is shown in Fig. [4j For /3 y = 1000, the turbulence 
survives at Rm = 3200 but dies at Rm = 1600, whereas for (3 y = 10000, the turbulence dies 
at Rm = 3200. It would seem that for each increase in (3 y by a factor of 10, the critical Rm 
value increases by roughly a factor of 2; it becomes easier to kill off the MRI with a lower 
resistivity as the background toroidal field is weakened. 



The time- and volume-averaged Elsasser numbers, calculated via equation (|8j), are 2.3 
and 137 for the (3 = 1000, Rm = 1600 and 3200 models respectively. The time average is 
from orbit 140 to the end of the calculation, and since the turbulence decays, the average A 
for Rm = 1600 agrees with that given by the initial toroidal field value. For the (3 = 10000, 
Rm = 3200 model, the Elsasser number is 0.45, again, calculated from either the initial 
toroidal field strength or from equation ^ after the turbulence has completely decayed. In 
both cases the initial poloidal field has f3 p = 1600, corresponding to an initial poloidal A on 
order unity for these Rm values. These results are consistent with the A values calculated 
in 



Simon & Hawley (2009); A < 100 leads to decay. 



The time- and volume-averaged Q y values are Q y = 34 and Q y = 33 for 32FTNum/31000 
and 32FTNum/310000, respectively. The time average is done from orbit 20 to 110. The 
toroidal field MRI appears to be quite well-resolved in both simulations. Averaging over 
the same time interval, we find that Q z = 7.3 for 32FTNum/31000 and Q z = 6.9 for 
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32FTNum/310000; the vertical field MRI is marginally resolved. Note that these Q values are 
calculated via the saturated state of these runs. Indeed, the initial Q values are substantially 
lower and either marginally or under-resolved. It is not clear which Q values are a better 
representation of how well-resolved the MRI actually is; the initial values correspond to the 
background field which ultimately drives the MRI, but other modes may become important 
in driving the MRI in the fully nonlinear state. In any case, despite the low initial Q values, 
sustained turbulence develops when explicit dissipation terms are not included, suggesting 
that at least some MRI modes are resolved. It is only when resistivity is turned on that 
decaying turbulence is observed. The initial low Q values may introduce uncertainty for the 
specific values of the critical resistivity, but the relationship between critical resistivity and 
initial field strength is likely to hold nevertheless. 



4. Vertically Stratified Simulations 



4.1. Baseline Simulations 



We now turn to a series of simulations that investigate the effects of physical dissipation 
in vertically stratified, isothermal disks. The various combinations of runs are summarized in 
Table [2] Some time- and volume-averaged values of various quantities are given in Table [3j 
The baseline calculations without physical dissipation are done at two resolutions: 32Num 
which uses 32 grid zones per H, and 64Num which uses 64. 32Num is run to a total 
integration time of 1058 orbits. It has sustained turbulence at a level of a = 0.028, calculated 



via equation (18), time-averaging from orbit 20 until the end of the calculation, and volume- 
averaging over all x and y and for \z\ < 2H. In Model 64Num, a = 0.022, where the time 
average is from orbit 20 until the end of the simulation at 159 orbits. 

One particularly useful diagnostic is the space-time diagram of horizontally averaged 
quantities. Figure [5] shows space-time plots of horizontally averaged B y and pressure- 
normalized total stress for a 100 orbit period in the 64Num simulation. From this, we 
see that B y changes sign periodically and rises above the equatorial plane. The behavior of 
By is roughly mirror symmetric around the equator and has a period of ~ 10 orbits at both 



resolutions. This has been seen in previous vertically stratified shearing boxes (e.g., Bran- 
denburg et al. 1995 Stone et al. 1996 Hi rose et al.|2006 ; Guan fc Gammie 2010 |Shi et aL 



2010 Gressel 2010 Davis et al. 2010), with a variety of initial fields, numerical resolutions 



and codes, and it has even been observed in global simulations (e.g., Fromang & Nelson 



2006 Dzyurkevich et al. 2010 O'Neill et al. 2010). It is apparently a generic property of 



MRI-driven turbulence in the presence of vertical gravity. 
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The white contours in the top panel of Fig. [5] indicate where the gas ft value switches 
from greater than to less than unity. For \z\ > 2 — 2.5H, ft < 1, except for some regions 
very near the vertical boundaries where ft > 1 as a result of the absence of magnetic field. 
The region around \z\ = 2H is also where the fluid becomes completely buoyantly unstable. 



Using the criterion of Newcomb (1961) as outlined in Guan & Gammie (2010), the gas is 
buoyantly stable if 



dp 

dz 



> 



P 2 9 



(19) 



where 7 = 1 because the gas is isothermal. According to this criterion applied to the 
simulation data, the fluid is buoyantly unstable for \z\ > 2H . Consistent with this, the 
space-time plot shows that the field rises faster for \z\ > 2H. For \z\ < 2H, there are 
regions of instability as well as stability, and it is within this region that the MRI is active as 
indicated by the presence of significant stress. Indeed, the total stress drops off rapidly near 
\z\ ~ 2H. It appears that the marginal buoyancy stability coupled with the MRI-induced 
turbulence leads to a slower rise of magnetic field until \z\ ~ 1.5-2H, where the gas then 
becomes completely buoyantly unstable. These results are consistent with the recent ZEUS 



calculations of Guan & Gammie (2010) with large radial extent as well as with Shi et al. 



(2010) using a version of ZEUS that includes radiation physics and total energy conservation. 



Compare, e.g., the top panel of Fig. 5 to Fig. 6 in Shi et al. (2010). 



Figure [6] shows the volume-averaged toroidal field within \z\ < 0.5H as a function of 
time. The RMS field fluctuations averaged within this region are roughly a factor of ~ 4 
times larger; the turbulent fluctuations dominate over the mean, but not significantly so. 
The temporal oscillations in the mean field are apparent. The right figure shows the tem- 
poral power spectrum, revealing the dominant ~ 10 orbit period. The oscillation amplitude 
appears to be modulated on longer timescales, ranging from tens to hundreds of orbits. 
Furthermore, the averaged radial field appears to exhibit the same 10 orbit cyclic behavior 
as the toroidal field, but with a slight temporal lag, as shown in Fig. [7| The behavior of 
the mean radial and toroidal fields resembles the a-Q dynamo model derived in |Guan 



Gammie (2010). Specifically, they write simplified evolution equations for volume-averaged 



field components, 



d(B y ) 
dt 



-qSl(B x ) 



\va\ 
2H 



(By) 



a>i 
2H 



(B x ), 



(20) 



d(B x ) 
dt 



\va\ 
2H 



(B x ) 



a 2 
2H 



(By). 



(21) 
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The first term on the right hand side of equation (20) is the background shear acting on the 
radial field. The second term is the loss of toroidal field due to buoyant rise, estimated to 
have a characteristic buoyant velocity equal to the toroidal field Alfven speed. The third term 



is the a-dynamo term coupling B x to the evolution of B y . Equation (21) is similar except 



there is no shear term and the toroidal and radial field components have been flipped with 



respect to equation (20). Also, in general, a\ ^ a 2 . Note that while we follow the notation 



and general approach of Guan & Gammie (2010) here, Gressel (2010) has also constructed 



an a-Vt model, which incorporates a helicity-based dynamo mechanism and reproduces the 
observed behavior of the horizontally averaged field components. 



Guan & Gammie (2010 ) numerically integrated this set of equations and found a solution 
that looks strikingly similar to the red and black curves in Fig. [7j The value of a\^ 2 sets 
the oscillation frequency in the mean field; ot\ = a 2 = —0.01QH reproduces the 10 orbit 



variability observed in stratified shearing boxes. Here we numerically integrate equation (20) 
using the simulation data for (B x ) (the red curve) and the initial condition for (B y ) taken 
from (By) at t = 0. We have set ot\ = so that the (B y ) evolution is controlled entirely by 
shear and buoyancy. The result is the blue curve in the figure. The agreement between the 
actual evolution of (B y ) and the a-Q dynamo model suggests that the evolution of the mean 
toroidal field within the mid-plane region is almost completely controlled by the shearing of 
radial field and the buoyant removal of the generated toroidal field. 

The remaining question then is, what creates the mean radial field? In the dynamo 
model, the a 2 term is essential, but what mechanism is responsible for a 2 ^ 0? The 
most likely candidate is MRI turbulence; turbulent fluctuations create EMFs that gener- 



ate poloidal field (e.g., Brandenburg et al. 1995 Davis et al. 2010 Gressel 2010), but the 
details of this are still not well understood. 

In Fig. [8j we plot the time- and horizontally-averaged vertical distributions for various 
quantities in 32Num. The time average runs from orbit 100 to the end of the calculation. 
The figure shows that the stress drops off rapidly near \z\ ~ 1.5H. The shape of the distribu- 
tion is generally the same for both Maxwell and Reynolds stresses, with the Maxwell stress 
always greater than the Reynolds stress by a factor that varies from 2.3 to 5.6 depending 
on z; this factor is ~ 4 when averaged over all z, in agreement with unstratified simulations 



(e.g., Hawley et al. 1995). The gas pressure maintains an approximately Gaussian distri- 
bution consistent with hydrostatic equilibrium, while the magnetic pressure is subthermal 
and relatively flat for all \z\ < 2H. The magnetic field becomes superthermal for \z\ > 2H, 
although its magnitude continues to decrease with height. These results are consistent with 
previous studies of isothermal disks (e.g 



Stone et al. 


1996 


Miller & Stone 


2000; 


Guan & 



Gammie 2010). Interestingly, the vertical structure of the turbulence is also consistent with 
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local simulations containing radiation physics (Hirose et al. 2006 Krolik et al. 2007 Hirose 



et al. 2009; Flaig et al. 2010) and global simulations (Hawley et al. 2001 Fromang & Nel- 



son 



2006), though we do not observe the double peak profile in the stress as seen in the 



simulations of Hirose et al. (2009) and Flaig et al. (2010). 



Figure [9] examines the three-dimensional structure of the magnetic field in the fully 
turbulent gas using streamline integration for 32Num at t = 100 orbits. The field strength 
is denoted by color rather than field line density. Within \z\ < 2H, the field is primarily 
toroidal but has a smaller scale, tangled structure in the x and z directions. Very near 
the vertical boundaries, however, the field appears to develop larger poloidal components. 
Utilizing snapshots taken throughout the evolution of 32Num, we find that this is typical of 
the saturated state, except for at t — 550 orbits, in which the field near the boundaries is 



primarily vertical. Figure 16 of Hirose et al. (2006) shows the field in a shearing box with a 
vertical domain twice as large as 32Num. In their simulation, the field at \z\ = 4H is mainly 
toroidal while the field structure at \z\ = 4H in 32Num resembles the field structure near 
\z\ = 8H in their simulation. This suggests that the vertical outflow boundary conditions are 
influencing the field structure very near the boundaries. Away from the vertical boundaries, 
however, the magnetic field in 32Num appears to have a very similar structure to that in 



Hirose et al. (2006). 



4.2. High and Low Turbulence States 

Having established the baseline simulations without physical dissipation, we now turn 
to the main focus of our present work: the effect of v and rj on vertically stratified MRI 
turbulence. 

The Maxwell and Reynolds stress time evolution for the simulations performed with 64 



zones per H is shown in the left panel of Fig. 10 Runs 64Rm3200Pm0.5, 64Rm800Pm0.5 



and 64Rm3200Pm4 have decreasing turbulence levels, while turbulence is sustained at a 
statistically constant value in the remaining simulations. Furthermore, 64Rm800Pm0.5 un- 
dergoes alternating periods of low and high stress, though the overall trend is downward 
with time. 



The right panel of Fig. [10] is the stress evolution for the equivalent simulations with 32 
zones per H, shown for a much longer time period than the 64 zone runs. There is con- 
siderable variability on long timescales. 32Rm3200Pm0.5 in particular exhibits alternating 
periods of low and high stress levels, occurring on ~ 100 — 200 orbit timescales. The more 
viscous and resistive run, 32Rm800Pm0.5, shows similar variability but on a much shorter 



-18- 



timescale of > 10 orbits. Turbulence in 32Rm3200Pm2 appears to have leveled off at a 
smaller value without any indication of regrowth to the higher state. Both 32Rm6400Pm4 
and 32Rm3200Pm4 remain at relatively high stress levels, which are very similar between 
the two runs, likely a result of having the same P m . 



Turning to the high and low states in 32Rm3200Pm0.5, the left panel of Fig. 11 shows 
the vertical profile for the horizontally averaged total stress during a high turbulence state 
(solid line, time averaged from 500 to 570 orbits) and a low turbulence state (dashed line, 
time averaged from 700 to 770 orbits). The high state has considerable stress within the 
mid-plane region (\z\ < 2H), whereas in the low state, the stress is peaked near \z\ ~ 2H. 
Note that there is still a nonzero stress in the mid-plane region even during the low state. 
The Reynolds stress is nonzero and relatively flat throughout the mid-plane, whereas the 
Maxwell stress drops off dramatically at the mid-plane. What Maxwell stress is present in 
the low states seems to derive from the presence of residual radial and toroidal field. (This 



effect was also noted by Turner et al. (2007) in the dead zone regions of their shearing box 
calculations). One characteristic of MRI turbulence is the ratio of the Maxwell stress to the 
magnetic pressure, which is typically on order 0.4 ( Hawley et al.fl995 Blackman et al.||2008 



and Table |3J). The vertical profile of this ratio is shown in the right panel of Fig. 11 This 
ratio is between 0.3 and 0.4 for \z\ < 1.5H in the high state and near \z\ ~ 2H in the low 
state; in the mid-plane region during the low state, the ratio is much smaller and approaches 
zero. 

These differences between the high and low states are seen in the other simulations that 
exhibit this variability (except for 32Rm800Pm0.5, in which the variability occurs on such 
a short timescale that temporal averaging becomes difficult). In summary, during the high 
state, the MRI is fully active, producing turbulent stresses for \z\ < 2H. In the low state, 
the turbulence in the mid-plane region has died down, leaving weaker Reynolds stresses. 
MRI-driven turbulence is present, however, in narrow regions near \z\ ~ 2H. This is similar 



to the behavior seen in Fleming & Stone (2003) where active layers above and below the 



mid-plane drive Reynolds stresses in the equatorial dead zone. 

Is the variability between high and low turbulence an artifact of using a relatively low 
resolution in these calculations? Several pieces of evidence suggest this is not the case. First 
of all, both resolutions for Re = 1600, P m = 0.5 show the same variable stress behavior on 
> 10 orbit timescales. Secondly, even at 32 zones per H, dissipation coefficients play a signif- 
icant role in determining the stress level, as shown above. The fact that the low resolution, 
unstratified simulations show constant saturation levels for sufficiently small 77 whereas verti- 
cally stratified simulations exhibit this variability for the same parameters suggests that the 
variability is a direct result of adding in vertical gravity. Thirdly, this variability was seen in 
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the simulations of Davis et al. (2010), which were run at a higher resolution of 64 zones per 



H. We note, however, that while Davis et al. (2010) used the same numerical algorithm as 



in this work, they used a different initial magnetic field configuration and vertical boundary 
conditions. 

Finally, we examine the Q values ( 10 ) using an Alfven speed defined by ^ and a volume 



average for all x and y and for \z\ < 0.5H. Figure 12 shows Q y and Q z as a function of time 



for 32Rm6400Pm4 and 64Rm6400pm4. The results suggest that the toroidal field MRI is 
quite well-resolved, but that the vertical field MRI may be only marginally resolved for the 
lower resolution simulation. The higher resolution Q values are roughly a factor of 2 larger 
than the lower resolution Q values, which is simply a result of the change in the grid zone 
size; the turbulent saturation level is roughly the same between the two resolutions. This 
result, coupled with the somewhat low Q z value, suggests that the vertical field MRI may not 
be playing a particularly significant role in setting the saturation level in these simulations. 

If this high-low variability is indeed a physical effect, what is its origin? First, consider 
the space-time diagram of the horizontally averaged B x and B y for several simulations. 



Figure 13 shows the first 200 orbits of 32Rm800Pm0.5. The turbulence level decreases 



dramatically from the beginning (see also Fig. 10). This is not too surprising considering 
that the same Rm, P m values cause rapid decay of the turbulence in the unstratified case; the 
resistivity is large enough to damp the MRI. The space-time plots show that after this decay, 
there is a residual magnetic field left within the mid-plane region of the disk. In particular, 
within \z\ < 0.5H, there is a net horizontally averaged B x < and B y < around t — 110 
orbits. The average B x within this region remains constant for awhile, and B y increases due 
to the shear of B x , eventually flipping to B y > 0. By t ~ 130 orbits, the turbulence has 
re-emerged, and the average B y rises to larger \z\. The resistivity then kills off the MRI 
again, leading to another period of B x shearing into B y before the next period of increased 
turbulence. 

Model 32Rm3200Pm0.5 also experiences alternating high and low turbulence states, 
but it is not immediately clear why the turbulence within the mid-plane region should 
decay since the unstratified run with the same dissipation terms has sustained turbulence. 



As noted earlier and in Simon & Hawley (2009), the critical Rm value below which the 
turbulence decays in unstratified shearing boxes is ~ 1000, but here Rm = 3200. The major 
difference between the stratified and unstratified simulations is that with vertical gravity, 
the net magnetic field within a localized region of the domain can change due to buoyancy, 
whereas with unstratified boxes the net toroidal field remains constant. Loss of flux can 
raise the critical Rm value (numerical resolution is also likely to have an effect); indeed, the 
unstratified simulations initialized with (3 = 1000 and (3 = 10000 fields have demonstrated 
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that the critical Rm depends upon the background field strength. The loss of net flux in a 
stratified box could similarly change the critical Rm value. 

In this simulation, the average toroidal field within the mid-plane region, (B y ), oscillates 
around zero with a period of 10 orbits. Thus, every 10 orbits or so, (B y ) is conceivably weak 
enough for resistivity to kill the turbulence. However, the turbulence remains for many of 
these 10 orbit periods. Furthermore, averaging B y within some vertical distance from the 
mid-plane erases information about the field structure there; e.g., (B y ) might be small but 
there could still be strong toroidal fields of opposite polarity close to z — 0. The point is, 
one cannot necessarily expect the turbulence to decay away strictly whenever (B y ) drops 
below a certain (small) value. 

The (B y ) oscillation amplitude appears to be modulated by a longer timescale, more 
on the order of ~ 100 orbits (see, e.g., Fig. [6]). This behavior is present in all simulations 
with and without physical dissipation. This is also roughly the same timescale over which 
the system switches between low and high states in the Rm = 3200 simulations. Comparing 
the time evolution of (B y ) for \z\ < 0.5H with the evolution of the total stress shows that 
the minima in the oscillation amplitude are generally correlated with the decay of mid-plane 
turbulence. One exception is near 200 orbits in 32Rm3200Pm0.5 in which (B y ) becomes 
rather small, but the turbulence remains active, though relatively weak compared to the 
fully active state. This correlation implies that if the mean toroidal field near the mid-plane 
remains sufficiently small for some time, resistivity can overwhelm the MRI and cause decay. 

Evidently, there exists a critical Rm below which the turbulence experiences long- 
timescale variability; this critical value is Rm < 6000. We carry out two more stratified 
simulations with Rm ~ 6000 but different P m values in order to further test this hypothesis. 
The first simulation is 32Rm6400Pm8; thus, Rm = 6400, and P m is relatively large. The 
mid-plane turbulence is sustained over a long time, nearly 330 orbits, without any sign of 
decay. The dissipation coefficients of the second simulation, 32Rm6250Pml, are chosen to 
match the relatively high Rm, low P m simulation that decays in the zero net flux shearing 



box (see Fromang et al. 2007; Simon & Hawley 2009). This simulation also remains in the 



high state for nearly 330 orbits. 

If the dissipation is large enough to cause MRI-driven turbulence to decay, what leads 



to its reactivation after > 100 orbits of time? Figure 14 shows the space-time plot of B x and 
B y for a 300 orbit period in 32Rm3200Pm0.5 during which the mid-plane turbulence dies 
out and eventually returns. For clarity, we also plot the volume-averaged horizontal field, 
(B x ) and (B y ), where the average is done for all x and y and for \z\ < 0.5H. During the low 
state, there is a small net radial field that remains in the mid-plane region and generates B y 
through shear. The sign of the net (B x ) reverses and with it the sign of d(B y )/dt. The last 
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such switch in the low state occurs around 750 orbits after which (B y ) continues to grow up 
until 790 orbits. At this point, the 10 orbit period oscillations resume. 



Figures 13 and 14 suggest that it is the growth of B y due to shear that periodically 
reactivates the MRI in the mid-plane. Indeed, during the low states, the poloidal fields are 
sufficiently weak that the most unstable wavelengths of the radial and vertical field MRI are 
very under-resolved; Q X <1 and Q z < 1 for both 32Rm800Pm0.5 and 32Rm3200Pm0.5. (Q 
was calculated as a function of time using equations ^ and (10).) Considering the toroidal 



field, however, we find that Q y is well above the marginal resolution limit when the mid-plane 
turbulence starts to decay. In the low states of 32Rm800Pm0.5, Q y ~ 10 — 20 occasionally 
dropping to Q y = 6. The typical Q y values for the low states of 32Rm3200Pm0.5 are similar 
but somewhat smaller. Of course when studying the behavior of the MRI in a situation where 
the Q values are small, one is in a regime where numerical resolution is likely to matter. 
The specific behavior of the disk in the low state might be different at higher resolution, but 
shearing of radial into toroidal field is itself not so dependent on resolution. Thus, Q may 
well play a role in determining when the MRI is reactivated, but not how it is reactivated. 

As a further test, we carry out two additional experiments. First, we take the state 
of the gas in 32Rm3200Pm2 at t — 150 orbits when the stress levels are decreasing and 
the average of B y within \z\ < 0.5H, (B y ) = 5.9 x 10~ 6 , is relatively small compared to 
the oscillation amplitude of (B y ) in the high state, which is ~ 5 x 10 -5 . We restart this 
simulation and add a net B y = 8.9 x 10~ 5 into the region \z\ < 0.5H, which corresponds to 
a toroidal /3 ~ 126 (using /3 defined with the initial mid-plane gas pressure P a ). This run is 
32Rm3200Pm2 _By+ (see Table [2]). Figure 15 shows the subsequent evolution of the stress 



along with the stress evolution of 32Rm3200Pm2. Not only does the mid-plane turbulence 
return, but the system undergoes episodic transitions between low and high states on ~ 100 
orbit time scales, as in 32Rm3200Pm0.5. 



In a second experiment, we initialize a stratified shearing box with all the same param- 
eters as in 32Num but with an initially very weak radial magnetic field. Specifically, for 
\z\ < 0.5H, B x = —<JlP l (5 X where /3 X = 10 6 . This field strength is very under-resolved; 



the Q x value is 0.2, and the radial field MRI will not be significant. Figure 16 shows the 
space-time diagrams of horizontally averaged B x and B y . The weak radial field is sheared 
into toroidal field that grows linearly in time until it reaches a sufficient strength to activate 
the MRI. The onset of MRI turbulence is rapid, and once the MRI sets in, the subsequent 
behavior is very similar to the other vertically stratified MRI simulations; there are rising 
magnetic field structures, dominated by the toroidal component, and the period of oscilla- 
tions in the mean field is ~ 10 orbits. Note that this experiment is not an exact imitation of 
the low state in our simulations; there is no MRI activity near \z\ ~ 2H initially. However, 
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this calculation demonstrates that shear amplification of a weak radial field can eventually 
lead to turbulence and move the system to the high state. 

All of the simulations in which turbulence in the mid-plane sets in after a period of 
quiescence show the presence of a net radial field in the mid-plane during the low state. 
32Rm3200Pm2, however, is the only simulation that does not show the re-emergence of the 
mid-plane MRI, despite over 1000 orbits of integration. An examination of the mid-plane 
region (up to \z\ < H) in the low state of this run shows that the residual radial field is 
weaker than in the low states of the other simulations. If this radial field would remain 
constant in time, the toroidal field would continually strengthen to the point of reactivating 
the MRI. However, (B x ) continues to change sign even in the absence of turbulence, though 
with a period of many hundreds of orbits. That is, (B x ) oscillates about zero but with a 
very small amplitude. (B y ) similarly oscillates around zero and never reaches a sufficient 
amplitude to reactivate the MRI, and the simulation remains in the low state. 

In fact, (B x ) oscillates about zero in all of our simulations, even in the low states (see 



e.g., Fig. 14). The various space-time diagrams suggest horizontal field migrates toward the 
mid-plane region from near \z\ ~ 2H where MRI activity is on-going. The oscillations in 
the mid-plane field components are a direct result of the oscillations previously observed 
in vertically stratified MRI turbulence. It is not entirely clear, however, what causes the 
field to migrate to the mid-plane. Does it diffuse there or is it carried downward by some 



sort of large scale flow? The shearing box calculations of Turner et al. (2007) and Turner & 



Sano (2008), which include Ohmic resistivity via a treatment of ionization chemistry, suggest 
that turbulent diffusion from the active layers is responsible for the field migration. While 
our simulations include a simpler prescription for resistivity, it is conceivable that a similar 
process is at work here. These issues will be addressed in future work. 

To summarize the results for transitions between high and low states, we find that for 
sufficiently small Rm, MRI turbulence within the mid-plane region can decay away, in part 
because of loss of net flux due to buoyancy. Residual net radial field remains and subse- 
quently creates toroidal field from shear. Once this toroidal field reaches a sufficiently large 
amplitude, the MRI is reactivated, and the turbulence is sustained for some duration until 
it decays again, repeating the pattern. In other words, the a-Q dynamo, seen in strati- 
fied simulations without resistivity, continues to operate in resistive simulations where the 
mid-plane MRI is suppressed. The dynamo is key to reactivation of MRI-driven turbulence. 

This behavior appears to be independent of the P m values we have used, except for near 
Rm = 3200; 32Rm3200Pm4 remains in the high state, whereas 32Rm3200Pm2, 32Rm3200Pm2_By+, 
and 32Rm3200Pm0.5 do not. However, in the higher resolution runs, 64Rm3200Pm4 and 
64Rm3200Pm0.5 show decay but 64Rm3200Pm2 appears to be sustained (though again, 
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these simulations were not integrated very far). While P m may play a role here, the line 
between sustained and highly variable stress levels is unlikely to be hard and fast, and many 
factors probably contribute to the nature of the turbulence near this Rm value. 



Comparing our results with those of Davis et al. (2010), we note that the largest Rm 



used in their simulations was Rm = 3200, consistent with the largest critical Rm observed 
in our simulations. Secondly, an examination of the space-time data from their simulation 
with Rm = 1600, P m = 2 (kindly provided by the authors) shows the same behavior as 
we have observed here; a net radial field remains within the mid-plane region after decay, 
shearing into toroidal field, from which the MRI is reactivated within this region. 



Interestingly, Turner et al. (2007) briefly notes that the toroidal field in the dead zone of 
some of their highly resistive simulations can occasionally become strong enough to reactivate 
turbulence there, after which the MRI is rapidly switched off again due to Ohmic dissipation. 
These simulations are consistent with our highest resistivity (Rm = 800) simulations, in 
which turbulence occurs in short outbursts rather than through a transition to a long-lived 
high state, as is the case at lower resistivities. 

Lastly, we examine the topology of the magnetic field in the shearing box during 
the low state. Figure 17 shows the equivalent information as Fig. [9j but for orbit 550 of 



32Rm3200Pm2_By+. For \z\ > 2H, the field remains mainly toroidal but with noticeable 
excursions from being purely toroidal, resembling the field structure in this region during the 
high state. Within 2H, the field is almost completely toroidal, and any small poloidal field 
present within this region is not visible in this image. We also examined the azimuthally 
averaged poloidal field structure in several snapshots of this run. We found that the struc- 
ture of the field was different, depending on which snapshot we examined. At some times, 
the poloidal field within 2H is almost completely radial, with very little vertical field. At 
other times, the vertical and radial fields are comparable in size such that the field takes on 
a more loop-like structure. 



4.3. The Prandtl Number Effect on Sustained Turbulence 

In this section, we investigate how P m affects angular momentum transport in simu- 



lations that do not exhibit the high-low variability. Figure [18] shows the volume-averaged 
stress levels for these simulations normalized by the gas pressure. The volume average is 
done for all x and y and for \z\ < 2H. The figure shows a general increase in the turbulence 
levels with P m , but there is also significant time variability in the stress. As a result, the 
curves overlap at certain times, much more so than for unstratified turbulence (see Fig. |2j). 
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Figure [19] displays the time-averaged a values as a function of P m for the unstratified 
(left panel) and stratified simulations (right panel). The time average for the unstratified 
simulations is the same as in Fig. [3j from orbit 120 to the end of the calculation, and for the 
stratified simulations, this average is done from orbit 150 until the end of the simulation. 
The error bars denote one standard deviation about the temporal average of the numerator 
in While there is a clear dependence of a on P m in the stratified simulations, it is not 
as steep as in the unstratified simulations. Taking a linear fit in log-log space, and assuming 
a oc P^, 5 = 0.54 for the unstratified runs (from § [3J, and 5 = 0.25 for the stratified 
calculations. In addition, the time variability is significantly larger in the stratified runs 
than in the unstratified case. 



Figure 20 is the vertical profile of the time- and horizontally-averaged total stress for 
the sustained turbulence runs. The time average was done from orbit 150 to the end of 
each simulation. Increasing P m increases the stress at nearly all z, and in all cases the stress 
drops off dramatically above \z\ ~ 1.5H, consistent with the baseline simulations. Run 
32Rm6400Pm8 appears to have a sharper peak in the stress profile near z = 0, compared 
to a flatter stress profile for \z\ < 1.5H in the other simulations. Using different temporal 
averaging windows for the averaged stress profiles produces the same general result; the 
stress increases with P m and 32Rm6400Pm8 has a sharper peak near the mid-plane. In some 
cases, however, the stress does not necessarily increase monotonically with P m at \z\ > 1.5H. 
Finally, we examined the vertical profile of the same quantities as in Fig. [8} We found that 
these profiles in the sustained turbulence simulations are all very similar to each other and 
to 32Num; the general vertical structure of the disk does not appear to be sensitive to P m . 

In summary, where the MRI is operating within a stratified disk, increasing P m leads to 
an increase in stress, consistent with the behavior seen in unstratified simulations. 



5. Summary and Discussion 



We have carried out a series of shearing box simulations with the Athena code to study 
MRI-driven turbulence with both vertical gravity and physical dissipation. Until the recent 
work of Davis et al. (2010), the role of physical dissipation in setting the level of angular 
momentum transport had been mainly studied in the context of unstratified simulations. 
As Davis et al. (2010) has shown, however, stratified simulations reveal new behaviors. 
Turbulence that decays in unstratified simulations is sustained in the presence of vertical 



lr rhe fluctuations in the volume-averaged pressure are very small and do not contribute significantly to 
the variability. 
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gravity and with large amplitude, ~ 100 orbit fluctuations in the stress in some cases. 



Our primary goal in this study has been a deeper investigation into the roles of viscosity 
and resistivity in MRI-driven turbulence using simulations that more accurately resemble as- 
trophysical disk systems. To this end, we have implemented more realistic, outflow boundary 



conditions in the vertical direction (Davis et al. (2010) use vertically periodic boundaries), 
and we have run the simulations for significantly longer integration times than is usual for 
shearing box calculations, in order to study long timescale effects. From these simulations 
we observe the following: 1) MRI turbulence is largely confined within 2 H of the equator. 
Above that height, the magnetic field becomes completely buoyantly unstable. 2) The mean 
field within the MRI-active portion of the disk evolves in a manner consistent with an a-VL 
dynamo. 3) Stratified disks show an increase in stress with increasing P m as previously seen 
in unstratified shearing boxes. 4) For Rm below some critical value, MRI turbulence in the 
mid-plane can be quenched putting the disk into a "low" state. 5) While in the low state, the 
continued action of the dynamo can raise the toroidal field to a level at which the mid-plane 
MRI is reactivated, switching the disk back to a "high" state. 

The first two of these results are consistent with previous stratified shearing box simula- 
tions and do not appear to be altered by the inclusion of physical dissipation. For \z\ < 2H 
in sustained turbulence simulations, the time- and horizontally-averaged turbulent energies 
and stresses are roughly constant with height, the magnetic field is only marginally stable to 
buoyancy, and /3 > 1. The predominantly toroidal field rises slowly away from the midplane. 
Above this height, the turbulence is significantly weaker, (3 < 1, and the field is completely 
buoyantly unstable, rising at a faster rate than for \z\ < 2H. These results, which are consis- 



tent with the ZEUS-based results of Guan & Gammie (2010) and Shi et al. (2010), suggest 



that there are two separate vertical regions in these disks. 

As a practical matter, this result suggests that to capture the behavior of the vertically 
stratified MRI, one need not go much beyond ±2H from the mid-plane. This is confirmed 
by a few additional simulations in which we extended the outer boundary to QH from the 
mid-plane instead of 4if . We found no difference in the vertical structure of the disk, the 
volume averaged stress levels, or the temporal variability of the system. 



Consistent with previous stratified shearing box simulations (e.g., Stone et al. 1996 



Hirose et al. 2006 Guan & Gammie 2010 Shi et al. 2010 Gressel 2010; Davis et al. 2010) 



and global simulations ( Fromang fc Nelson|2006 Dzyurkevich et al.|2010 O'Neill et al.|2010 ), 
we found a strong 10 orbital period variability in the mean magnetic field within the mid- 
plane region. This variability is characterized by the buoyant rise of predominantly toroidal 
field from the mid-plane region, after which a field of opposite sign takes its place at the 
mid-plane. The evolution of this toroidal field can be well-modeled by a simple a-Q mean 
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field dynamo where the shear of radial field into toroidal field and the subsequent buoyant 
removal of toroidal field from the mid-plane dominate the evolution. 

In addition to the dynamo period of roughly 10 orbits, we also observe a longer ~ 100 
orbit timescale variability. This behavior has not previously been reported, and while we do 
not know its origin, there is evidence that it plays a role in another novel effect resulting from 
vertical gravity: high and low states. For disks that are sufficiently resistive, i.e., have low 
enough Rm, MRI turbulence decays near the mid-plane leading to a quiescent period that 
is eventually followed by regrowth of the MRI in this region. The period of this variability 
is on a similar timescale to the long-period variability seen in all the stratified simulations. 
The transition to a low state occurs when the mean toroidal field in the mid-plane region 
is reduced below some critical value, presumably by buoyancy. Resistivity dominates over 
the MRI and the turbulence decays. An extremely weak mean radial field remains however. 
Although this mean radial field itself can vary during the low state, it creates toroidal field 
through shear. Shear amplification is largely unaffected by resistivity. Once the toroidal 
field reaches a particular strength, the mid-plane MRI is re-energized, and the disk becomes 
fully turbulent again, returning to the high state. 

This behavior is not particularly sensitive to P m and appears to be the same as that 
reported by Davis et al. (2010). Thus, it is not likely to be related to the dynamo issue 
of P m ~ 1 in zero net magnetic flux shearing boxes, which Davis et al. (2010) specifically 



investigated. Instead, it appears to be a robust behavior that emerges whenever the disk is 
sufficiently resistive. The critical Rm below which the high-low variability exists is 3200 < 
Rm c < 6000. If Rm > Rm c , the turbulence remains sustained for the dissipation parameters 
explored here, and averaged stress increases with P m for all z, though with a less steep 
dependence of a on P m compared to unstratified simulations. 

Because we are dealing with the decay of turbulence and the behaviors of weak magnetic 
fields, the values of properties such as Rm c may well be resolution dependent. However, 
given the long evolution times probed here, it was necessary to limit the resolution to 32 
grid zones per H. Before running the vertically stratified simulations, we executed a series of 
unstratified calculations at this resolution to test whether or not the MRI and the effects of 
physical dissipation would likely be resolved. We found that the effect of physical dissipation 
on MRI saturation is reasonably well converged at 32 zones per H: a increases with P m at 
this resolution, although with a steeper dependence than at higher resolution, and a is 
roughly constant with resolution above 32 zones per H for a given set of Re, Rm, and P m 
values. This in itself is an interesting result since it shows that the higher resolutions used in 



Fromang et al. (2007), Lesur & Longaretti (2007), and Simon & Hawley (2009) may not be 



necessary to capture the general effects of physical dissipation. It is also in agreement with 
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the recent vertically stratified shearing boxes of Flaig et al. (2010) that included radiation 
physics; they found that only ~ 30 zones per H in the vertical dimension are required to 
gain a reasonable representation of the turbulent saturation level. Thus, while the exact 
value of Rm c may change at higher resolution, the general effects of physical dissipation on 
the vertically stratified MRI are likely captured at the resolutions we used. 

What do these results imply for the MRI and its application to astrophysical disks? 
First, it had been previously appreciated that sufficient resistivity could cause a transition 
from a turbulent ("high") state to a non-turbulent ("low") state. Here we have found that 
a process resembling an a-Q dynamo can accomplish the reverse and re-establish the high 
state. The transition occurs on longer timescales than the ~ 10 orbit associated with the 
dynamo. This temporal variability could have potential applications for several types of 
accretion disks. Protoplanetary disks have large regions of low ionization gas, including a 



dead zone layer (Gammie 1996) where resistivity is too high to sustain the MRI. Dwarf 



nova disks also contain regions of partial ionization, and it is intriguing that the Rm values 
in these systems are on the same order as the critical Rm ~ 10 3 for the decay /regrowth 



behavior (Gammie & Menou 1998). Even some regions of AGN disks may have moderately 



high resistivity, though typical Rm values are probably larger than those in dwarf nova 



systems because of the larger disk scale height (Menou & Quataert 2001) 



It is tempting to associate the peaks and dips of turbulent activity in our simulations 
with the outbursts and variability observed in these systems. Of course, there remains much 
work before such a connection can be substantiated. In particular, more realistic simulations 
will have r] (and v) that depend on temperature and density, rather than being constant 
throughout the disk. Furthermore, the influence of other non-ideal MHD effects on the MRI 
needs more work. The Hall effect is often times just as important as Ohmic resistivity in 
astrophysical environments ( Wardle||l999 ; Balbus & Terquem 2001 Balbus|[2003 ), and while 
simulations including both Hall and Ohmic terms have been carried out by Sano fc Stone 



(2002a) and Sano & Stone (2002b), there remains more parameter space to explore and 



physics to include. Lastly, we note that if rj is so large that no MRI modes are present, there 
would be no temporal variability and the turbulence would be completely quenched. 

In disks that have Rm above the critical value, the MRI-driven turbulence operates 
continuously. Nevertheless, as previously seen in unstratified simulations, dissipation terms 
can have an impact: the stress level increases with increasing P m . This may be relevant to 
hot, fully ionized disk gas such as in X-ray binaries and some regions of AGN disks where P m 



has a strong temperature dependence (Balbus & Henri 2008). While our simulations show 



that angular momentum transport increases with P m , the Re and Rm values of such disks 
are significantly larger than the values probed here. Whether or not the P m effect continues 
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into the large Re/Rm regime remains very much an open area of research. 



Finally, one particular field geometry that has not been explored here or in most verti- 
cally stratified local simulations is that of a net vertical field. These simulations are quite 



challenging; the channel mode dominates the solution (Miller & Stone 2000 Latter et al. 



2010), leading to very strongly magnetized regions of the disk that can often times cause the 



numerical integration techniques to fail (but see Suzuki & Inutsuka (2009), in which a stable 
evolution was produced). 

In summary, we have explored the spatial and temporal behavior of the MRI in the 
presence of both vertical gravity and physical dissipation. We find that for moderately 
resistive simulations, the disk cycles between states of low and high turbulent stresses, and 
that orbital shear of radial field into toroidal field is essential to both this behavior as well as 
the temporal variability of sustained turbulence. When sustained, the stress increases with 
P m , in agreement with unstratified simulations. Our calculations are an important stepping 
stone towards more realistic simulations that include temperature-dependent v and rj. 



We thank Xiaoyue Guan, Juilan Krolik, Mitch Begelman, Phil Armitage, and Rosalba 
Perna for useful discussions and suggestions regarding this work, and we are very grateful to 
Shane Davis, Jim Stone and Martin Pessah for providing us with some of the data from their 
paper. This material was supported by NASA Headquarters under the NASA Earth and 
Space Science Fellowship Program - Grant NNX08AX06H; NASA grants NNX09AG02G, 
NNX09AB90G, and NNX09AD14G; by the NSF under grants AST-0807471 and AST- 
0908869; and by a Virginia Space Grant Consortium (VSGC) fellowship. The simulations 
were run on the TACC Sun Constellation Linux Cluster, Ranger, supported by the National 
Science Foundation. 



A. Numerical Methodology 



A.l. Integration 

The algorithms for solving the shearing box equations Q-(|3]) and the specific imple 
mentation of those algorithms within the Athena code are described in Stone Gardiner 



(2010); here we provide a brief summary of the numerical methods detailed in that paper. 



The left hand sides of the equations are solved via the standard Athena CTU flux- 



conservative algorithm (Stone et al. 2008). The gravitational and Coriolis source terms in 



the momentum equation are evolved via a combination of an unsplit method consistent 
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with CTU and Crank-Nicholson differencing constructed to conserve energy exactly within 



epicyclic motion (Gardiner & Stone 2005; Stone & Gardiner 2010). 



For our simulations in which the radial size of the shearing box is a scale height, H, the 
velocity is initialized with a background shear flow given by 



-qQx. 



(Al) 



For sufficiently large x domains, this velocity can become supersonic (orbital velocities are, 
in general, supersonic in accretion disks). For large x domains, the Courant limit on the 
timestep can become significant. Furthermore, the background shear flow can lead to a 
systematic change in truncation error with radial position in the box, which in turn introduces 



purely numerical features in the radial density and stress profiles (Johnson et al. 2008). 



Hence, for large radial domain simulations we utilize an orbital advection scheme, which 
subtracts off this background shear flow and evolves it separately from the fluctuations in 



the fluid quantities (Masset 2000; Johnson et al. 2008 Stone & Gardiner 2010). 



A. 2. Boundary Conditions 



In the shearing box, the y boundary condition is periodic. The x boundary condition is 
shearing-periodic (Hawley et al. 1995 Stone fc Gardiner]|2010 ): quantities are reconstructed 
in the ghost zones from appropriate zones in the physical domain that have been shifted along 
y to account for the relative shear from one side of the box to the other. In Athena, this 
reconstruction step is performed on the fluid fluxes, and then the ghost zone fluid variables are 
updated via these reconstructed fluxes. The order of this reconstruction matches the spatial 
reconstruction in the physical grid, e.g., 3rd order reconstruction of the ghost zone fluxes is 
done when the PPM spatial reconstruction is employed. Furthermore, the y momentum is 
adjusted to account for the shear across the x boundaries as fluid moves out one boundary 
and enters at the other. The y component of the electromotive force (EMF) is reconstructed 



at the radial boundaries to ensure precise conservation of net vertical magnetic flux (Stone 



& Gardiner 2010). 



In unstratified shearing box simulations the z boundary is periodic. This same boundary 



condition has been employed in stratified boxes as well (e.g. Davis et al. 2010). For our 
simulations, we use an outflow boundary condition instead. The density p is extrapolated into 
the ghost zones based upon an isothermal, hydrostatic equilibrium, using the last physical 
zone as a reference. Therefore, for the upper vertical boundary, the p value in grid cell k is 
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p(k) = p(ke)e W (- Z(fc)2 ^ (fce)2 ) , (A2) 

where ke refers to the last physical zone at the upper boundary. An equivalent expres- 
sion holds for the lower vertical boundary. This extrapolation provides hydrostatic support 
against the opposing vertical gravitational forces, which are also applied in the ghost zones. 
All velocity components, B x , and B y are copied into the ghost zones from the last physical 
zone assuming a zero slope extrapolation. If the sign of v z in the last physical zone is such 
that an inward flow into the grid is present, v z is set to zero in the ghost zones. Finally, the 
ghost zone values of B z are calculated from B x and B y to ensure that V • B = within the 
ghost zones. 



A. 3. Physical Dissipation 

Viscosity and resistivity are added via operator splitting; the fluid variables updated 



from the main integration (§ A.l) are used to calculate the viscous and resistive terms. The 
viscosity term is calculated via the divergence of the viscous stress tensor, equation Q, and 
the resistive term is included as an additional EMF within the induction equation (|3]). This 
formulation allows us to discretize the viscous and resistive terms in a flux-conservative and 
constrained-transport manner, consistent with the Athena algorithm. Note that the resistive 
contribution to the y EMF must also be reconstructed at the shearing-periodic boundaries 
in order to preserve B z precisely. 

The addition of viscosity and resistivity places an additional constraint on the timestep, 

( A 2 A 2 \ , s 

At = MIN lAt CTV ,C —,C —j , (A3) 



where AtcTU is the timestep limit from the main integration algorithm (see Stone et al. 2008 ), 
C Q = 0.4 is the CFL number, and A is the minimum grid spacing, A = MIN(Ax, Ay, Az). 
Note that in most of our simulations v and rj are sufficiently small that they do not limit 
the timestep. 
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Table 2. Vertically Stratified Simulations 



Label 


Re 


Rm 




Resolution 
(zones per H) 


Integration time 
istop - tstart (orbits) 


Description 


32Num 








32 


1058 


num. dissipation 


32Rm800Pm0.5 


1600 


800 


0.5 


32 


263 




32Rm3200Pm0.5 


6400 


3200 


0.5 


32 


863 




32Rm3200Pm2 


1600 


3200 


2 


32 


1082 




32Rm3200Pm4 


800 


3200 


4 


32 


487 




32Rm6250Pml 


6250 


6250 


1 


32 


337 




32Rm6400Pm8 


800 


6400 


8 


32 


325 




32Rm6400Pm4 


1600 


6400 


4 


32 


488 




32Rm3200Pm2_By+ 


1600 


3200 


2 


32 


584 


By added at 150 orbits 


32ShearBx 








32 


45 


net B x within midplane 


64Num 








64 


159 


num. dissipation 
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1600 
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0.5 


64 


108 
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64 
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64 


80 
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Fig. 1. — Time- and volume- averaged stress parameter a as a function of grid zones per H in 
the unstratified SZ simulations; a = ((pv x 5v y — B x B y )) / ((P)), where the average is calculated 
over the entire simulation domain and from 20 orbits to the end of the simulation. Only 
simulations with sustained turbulence are plotted. The squares are runs with Rm = 12800, 
P m = 16; asterisks are Rm = 12500, P m = 4; and triangles are Rm = 25600, P m = 2. By 32 
zones per H, the a values appear to be relatively close to the higher resolution values. 
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Fig. 2. — Volume-averaged total stress normalized by the the volume-averaged gas pressure 
as a function of time in orbits in the unstratified, toroidal field FT simulations. The black 
line corresponds to Rm = 3200 and P m = 4, magenta is Rm = 6400 and P m = 8, green is 
Rm = 800 and P m = 0.5, light blue is Rm = 3200 and P m = 2, dark blue is Rm = 6400 
and P m = 4, and red is Rm = 3200 and P m = 0.5. There is a clear increase in stress with 
increasing P m . For sufficiently low Rm, the turbulence decays. 
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Fig. 3. — Time- and volume-averaged stress parameter a as a function of P m in the unstrat- 
ified FT simulations (asterisks) and the higher resolution, net toroidal field simulations of 



Simon & Hawley (2009) (diamonds). In the FT simulations, a = ((pv x 5v y — B x B y ))/((P)), 



whereas in the higher resolution simulations, a = ({pv x 5v y — B x B y )) / P ; see Simon & Hawley 
(2009). These definitions are roughly equivalent since (P) ~ P a . For the FT simulations 



the average is calculated over the entire simulation domain and from 120 orbits to the end 
of the simulation. Only simulations with sustained turbulence are plotted. The dashed lines 
are linear fits to the data in log-log space. Both resolutions show a clear P m dependence, 
but this dependence is steeper at the lower resolution. 
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Fig. 4.— Volume-averaged total stress normalized by the volume-averaged gas pressure 
as a function of time in orbits for a series of unstratified shearing box simulations. The 
volume- average is done over the entire simulation domain. In each plot, the black line is 
from a simulation with only numerical dissipation. The colored lines are simulations with 
physical dissipation, initiated from the numerical dissipation run at orbit 100. The red lines 
correspond to Rm = 3200, P m = 2 and the blue line is Rm = 1600, P m = 1. The top panel 
is initiated with a background toroidal field characterized by /3 = 1000, and the bottom 
panel has = 10000. The weaker toroidal field appears to be killed off at a lower resistivity 
compared to the stronger background field. 
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Fig. 5. — Space-time diagram of the horizontally averaged B y (top panel) and total stress 
normalized by the volume-averaged gas pressure (bottom panel). The volume-average is 
done for all x and y and for \z\ < 2H. The white contours on the top panel denote where 
goes from greater to less than unity. The horizontally averaged B y appears to rise vertically 
into the upper z layers, being replaced in the mid-plane region by B y of the opposite sign. 
The rise speed of the field increases after \z\ ~ 2H is reached. The sign flipping in B y has a 
period of ~ 10 orbits. 
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Fig. 6. — Left: Time evolution of volume-averaged B y in 32Num. The volume average is 
done for all x and y and for \z\ < 0.5H. The dashed line corresponds to (B Xjy ) = 0. Right: 
Temporal power spectrum of (B y ) from the left plot, calculated from orbit 50 to 1050. The 10 
orbit period oscillations in (B y ) are immediately apparent in both plots, particularly as the 
peak in the power spectrum. The 10 orbit oscillations are modulated on longer timescales, 
ranging from tens to hundreds of orbits. 




Fig. 7. — Time evolution of volume-averaged field components for part of 32Num. Red is 
(B x ), black is (B y ), and blue is (B y ) as calculated from (B x ) using a simple a-Q dynamo 
model discussed in the text. The volume average is done for all x and y and for \z\ < 0.5H. 
The dashed line corresponds to (B x>y ) = 0. (B x ) has been multiplied by a factor of 5 relative 
to (By) to make a more direct comparison possible. The variations in (B x ) are accompanied 
by variations in (B y ), which are offset in time, and the dynamo model shows that the 
evolution of (B y ) is controlled by shear of radial field and buoyant removal of toroidal field. 
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Fig. 8.— Time- and horizontally-averaged vertical distributions for various quantities 
in 32Num. Upper left: gas density; upper right: Maxwell (solid) and Reynolds (dashed) 
stresses; lower left: gas pressure (solid), magnetic energy (dashed), and kinetic energy (dot- 
ted); lower right: gas (3 defined as the time- and horizontally-averaged gas pressure divided 
by the time- and horizontally-averaged magnetic energy density. The time average is done 
from orbit 100 to the end of the simulation. The stress and magnetic energy are relatively 
flat for \z\ < 1.5H but drop off substantially for larger \z\. Outside of \z\ ~ 2H, the magnetic 
energy dominates over gas pressure. 
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Fig. 9. — Magnetic field structure at t — 100 orbits in 32Num, produced via a stream line 
integration. The field strength (in code units) is displayed via color and not the density 
of the field lines. The magnetic field has a primarily toroidal structure but has a smaller, 
tangled structure in the x and z directions. 
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Fig. 10. — Volume-averaged total stress normalized by the volume-averaged gas pressure as 
a function of time in orbits in the vertically stratified simulations at 64 (left) and 32 (right) 
grid zones per H. The volume-average is done for all x and y and for \z\ < 2H. The black 
line corresponds to Rm = 3200 and P m = 4, green is Rm = 800 and P m = 0.5, light blue 
is Rm = 3200 and P m = 2, dark blue is Rm = 6400 and P m = 4, and red is Rm = 3200 
and P m = 0.5. Some of the simulations appear to undergo periods of low stress followed by 
higher stress, occurring on very long timescales of ~ 100 orbits in some cases. 
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Fig. 11.— Left: Time- and horizontally-averaged total stress as a function of z for 
32Rm3200Pm0.5. The stress is normalized by the time- and volume-averaged gas pres- 
sure, where the volume average is done for all x and y and for \z\ < 2H. Right: Ratio of 
time- and horizontally-averaged Maxwell stress to time- and horizontally-averaged magnetic 
energy in 32Rm3200Pm0.5. In both plots, the solid line corresponds to a time average from 
500 to 570 orbits, which is a state of high turbulence, and the dashed line is a time average 
from 700 to 770 orbits, a low turbulence state. During the high state, the stress is relatively 
flat for \z\ < 2H. In the low state, the stress is smaller within the mid-plane region and 
peaks near \z\ ~ 2H. 
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Fig. 12. — Quantitative measurement of how well resolved the MRI is in the toroidal (Q y ; left 
plot) and vertical (Q z \ right plot) directions as a function of time. The solid line corresponds 
to run 32Rm6400Pm4 and the dashed line is 64Rm6400pm4. The dotted horizontal line 
corresponds to Q — 6, below which the MRI is considered to be under-resolved (Sano et al. 



2004). Qi is calculated using the volume average of the magnetic energy and gas density 

z\ < 0.5H (see text). The toroidal field MRI at both resolutions as 



for all x and y and for 

well as the vertical field MRI at the higher resolution appear to be reasonably well-resolved. 
However, at the lower resolution, the vertical field MRI is only marginally resolved. 
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Fig. 13. — Space-time diagram of the horizontally averaged B x (top panel) and B y (bottom 
panel) for the first 200 orbits of 32Rm800Pm0.5. The turbulence initially decays, leaving 
a net B x within the mid-plane region, which shears into toroidal field. This appears to 
eventually reenergize the MRI, but the large resistivity quickly quenches the turbulence 
again. 
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Fig. 14. — Space-time diagram of the horizontally averaged B x (top left) and B y (bottom 
left) for a 300 orbit period in 32Rm3200Pm0.5, and the average of B x (top right) and B y 
(bottom right) over all x and y and for \z\ < 0.5H as a function of time in orbits for the same 
300 orbit period. During the period of no MRI turbulence, a net radial field still exists within 
the mid-plane region. This field appears to flip signs occasionally, leading to corresponding 
flips in B y due to shear. 
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Fig. 15. — Volume-averaged total stress normalized by the volume-averaged gas pressure as 
a function of time in orbits in 32Rm3200Pm2 (solid line) and 32Rm3200Pm2_By+ (dashed 
line). The volume average is done over all x and y and for \z\ < 2H. The run in which a net 
B y is added into the mid-plane region (dashed line) remains in the high state at first and 
then exhibits variability between low and high states. 




Fig. 16. — Space-time diagram of the horizontally averaged B x (top panel) and B y (bottom 
panel) for 32ShearBx. The uniform radial field that is present initially shears into toroidal 
field, which eventually becomes strong enough to launch the MRI. 
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Fig. 17. — Magnetic field structure at t = 550 orbits in 32Rm3200Pm2_By+, produced via 
a stream line integration. The field strength (in code units) is displayed via color and not 
the density of the field lines. The color scale is the same as that in Fig. [9] for comparison. 
This snapshot corresponds to the low state. The magnetic field is almost completely toroidal 
near the mid-plane but has a tangled structure in the upper \z\ regions, more reminiscent of 
the active state. 
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Fig. 18. — Volume-averaged total stress normalized by the volume-averaged gas pressure 
as a function of time in orbits in the lower resolution, vertically stratified simulations where 
turbulence remains sustained. The volume-average is done for all x and y and for \z\ < 2H. 
The left plot is the first 120 orbits of the evolution, whereas the right plot is 350 orbits of the 
evolution. The black line corresponds to Rm = 3200 and P m = 4, dark blue is Rm = 6400 
and P m = 4, magenta is Rm = 6400 and P m = 8, and brown is Rm = 6250 and P m = 1. 
The vertical axis has been chosen to match that of Fig. [2] for comparison. While the stress 
levels generally increase with P m , there is significant overlap between the different curves at 
different times. 
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Fig. 19.— Time- and volume-averaged stress parameter a as a function of P m in 
the unstratified FT simulations (left plot) and the stratified simulations (right plot); 
a = {{pv x Sv y — B x B y ))/((P)). The average is calculated over the entire domain (all x 
and y and for \z\ < 2H) and from 120 (150) orbits to the end of the simulation for the un- 
stratified (stratified) runs. The dashed lines are linear fits to the data in log-log space, and 
the error bars denote one standard deviation about the temporal average of the numerator in 
a. Both cases show a clear P m dependence. However, in the stratified runs, this dependence 
is less steep, and there is considerable temporal variability. 
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Fig. 20.— Time- and horizontally-averaged total stress as a function of z on a linear 
(left) and logarithmic (right) vertical scale. The stress is normalized by the time- and 
volume-averaged gas pressure, where the volume average is done for all x and y and for 
\z\ < 2H. The time average is done from orbit 150 until the end of each simulation. The 
solid line corresponds to 32Rm6400Pm8, the dashed line is 32Rm3200Pm4, the dotted line 
is 32Rm6400Pm4, and the dot-dashed line is 32Rm6250Pml. The stress appears to increase 
with P m for nearly all z, and for all P m , there is a sharp decrease in the stress for \z\ > 1.5H. 



